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Abstract 

We consider a supply chain network model consisting of buffering queues and proces- 
. sors. Such model was studied in |21] using a system of partial differential equations and 

ordinary differential equations. In this article, we propose a variational approach as an 
alternative to formulate this problem, and derive an algorithm for numerical computation. 
■ Based on this numerical method, a class of network control and optimization problems 

can be formulated as mixed integer programs (MIP). We conduct theoretical analysis of 
our algorithm and compare it with existing ones |20[ 121) , demonstrate its modeling and 
computational advantages which are then verified numerically. 

1 Introduction 

\Q ■ 1.1 Modeling overview 

Manufacturing systems can be described by using a number of mathematical approaches, which 
provide powerful tools to study and analyze the behavior of production systems under specific 
conditions. Among theses mathematical approaches, we distinguish between static/stationary 
models and dynamic models which are time-dependent, the latter will be the main focus of 
this paper. 

^ ■ The dynamic models describe and predict time evolution of system states by introducing 

dynamics to different production steps. These models can be further categorized as dis- 
crete event simulation (DES) and differential equation-based continuous models. The DES 
[5] is an exact and computationally intensive method, in which the evolution of the system is 
viewed as a sequence of significant changes in time, called events, for each part separately. A 
cost-effective alternative to discrete event models is fluid-like network modeling using partial 
differential equations, these continuous models were developed originally for traffic flow and 
have been studied extensively, for example, in [TH Q2J ISSJ ESI [281 E2]- Recently the ideas 
of continuous traffic models have been extended to supply chains, we refer the readers to 
[B El O HI [ini IS] for an overview. In particular, based on a simple rule for releasing parts, [2] 
derived a hyperbolic conservation law for the part density and flux in the supply chain: 

d t p{t, x) + d x m.in{V(x)p(t, x), fi(x)} = (1.1) 
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where p(t, x) denotes density of parts, V(x), fi(x) denote the processing speed and flow capac- 
ity, respectively. The above PDE will be asymptotically valid in regimes with a large number 
of parts and processors in the supply chain. It should be noted that the solution to the con- 
servation law (jl.ip can only be considered in the distributional sense, due to the discontinuous 
dependence of the flux function on x. This is easily seen from an example involving bottleneck: 
consider the flow capacity fi(-) at point a, and assume n{a~) > fi(a + ), if the maximum flow 
n(a~) is attained, a ^-distribution in the density p(t, •) will occur at x = a corresponding to 
a bottleneck. 

To overcome this difficulty, |21| uses a separate modeling of queues in front of each pro- 
cessor, avoiding ^-distribution in the density. This finally leads to a system of conservation 
law coupled with ordinary differential equations. This supply chain model can be easily ex- 
tended to incorporate general network topology and real-world production features such as 
multi-commodity or due-date production, as well as a wide range of mathematical problems 
such as network control and optimization |16j . 

1.2 Solution methods 

In [161120 1 [2T | I22|. the numerical methods for solving the system is an upwind finite difference 
discretization for the conservation law and a forward Euler scheme for the ordinary differential 
equation. To guarantee the stability of explicit discretization scheme, the CFL condition must 



The solution method described above encounters several numerical difficulties. First, the ODE 
in the system is discontinuous in its state variable, and several modifications were proposed 
to ensure the numerical performance [U[2D]. Second, despite these remedies, the system could 
still yield non-physical solutions in certain cases, see Section 15.51 for an example. Third, the 
CFL condition for the PDE and stiffness condition for the modified ODE [4J could potentially 
lead to large-size systems, due to constraints on time step size. 

In view of the above difficulties, we propose in this article a reformulation of the system 
using the integrals of flow (sometimes known as N-curves |28j) and Hamilton- Jacobi equations. 
This alternative not only eliminates the numerical problems above, but also leads to a more 
efficient algorithm and better solution precision. The corresponding solution algorithm is 
known as the variational method or Lax formula-based method and was developed in [61 Q31 
El [EJ [25] , it was originally proposed as a semi-analytic solution to the scalar conservation law 
and Hamilton- Jacobi equation [T71 [25], the application to fluid based models such as traffic 
flows are recently investigated in [H El El El [14] . Using the variational approach, the viscosity 
solution of the Hamilton-Jacobi equation is formulated as an optimization problem, which, 
depending on the PDE, may be simplified or explicitly instantiated. This provides us with 
a powerful alternative to existing finite difference-based numerical schemes. To list a few 
benefits of using the variational approach: 

(1) The Lax formula yields exact value of the solution to the PDE. 

(2) There is no need to discretize spacial domain and the values of the N-curves are given 
explicitly, this grid-free method allows for a faster computation and less memory consumption. 

(3) The algorithm does not impose any constraints on the time step or uniformity of the time 
grid, unlike finite difference schemes. 

(4) The resulting system is independent of spacial variable, this is arguably more general than 
the PDE-based system [20} [21] in terms of modeling assumption and solution method. 



hold [20] 
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An interesting application of the variational method proposed here is the formulation of 
network optimization problems as mixed integer programming (MIP). The MIP is one of the 
main formulations of supply chain optimization in the economic literature, see \30\ I34| . its 
connection with continuous network models brings new aspect to this type of problem and 
has been investigated in the context of traffic flows and supply chains |19[ [201 135] . This article 
reveals a new relation between MIP and fluid-based models from the viewpoint of variational 
principle. 

It is natural to compare the MIP derived in this article to the existing ones |19[ [20] . The 
latter is based on a finite difference discretization of the PDEs and ODEs. For the same reason 
mentioned before, the variational method-based MIP allows more efficient computation and 
yields better solution precision, compared to finite difference-based MIP. To obtain solutions 
with the same quality, the MIP we're presenting here will be significantly smaller in size, this 
will be established both theoretically and numerically in this paper. 

The rest of the article is organized as follows, in Section El we review the supply chain 
network model in [21] and discuss briefly the governing PDEs and ODEs. In Section [31 a 
variational approach is proposed to reformulate the system without invoking any PDEs or 
ODEs, the solutions are derived in both continuous and discrete time. Section H] considers the 
network optimization problem based on variational approach and derives a mixed integer pro- 
gramming, which is then compared to the MIP in [20] . Finally, several numerical experiments 
are presented in Section El 

2 Supply chain network model 

We begin our discussion with the network model presented in [21], it consists of separate 
modeling of queues (ODE) and processors (PDE). A precise description of the network is 
made via the notion of directed graph with edges (or arcs) and vertices (or nodes). 

Definition 2.1. (Supply chain network definition) 

1. A supply chain network is a directed graph (A, V) where each edge e G A corresponds 
to an individual processor, each vertex v £ V represents the respective queue. 

2. Each processor e G A is mapped onto an interval [a e , b e ], with b e — a e = L e being the 
length of the processor. 

3. Each processor possesses a queue, which is located at the vertex in front of the processor. 
4 ■ The flow capacity pf , processing speed V e and throughput time T e = L e /V e of each 

processor e £ A are constants. 

For each e G A, p e (t, x) denotes the density of goods at time t and location x G [a e , b e ], 
q e (t) denotes the size of queue in front of this processor. Assume a supply of goods at rate u e {t) 
enters the buffer queue, if any, before they're released to the processor (u e (t) is determined by 
the topology and flow profiles of the network, we'll postpone the details to Section f3.3j) . The 
dynamics in the processor and queue are governed by the following advection equation (|2.3p 



and ODE Q -flZSp. 

d t p e (t, x) + V e 8 x p e (t, x) = 0, p e (0, x) = p e (x), {t, x) G [0, +oo) x [a e , b% (2.3) 





(2.4) 




q e (t) = 0, 
q e {t) > 0. 



(2.5) 
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(|2.5|) represents the rate at which goods in the queue are released into the processor, this is 
sometimes called the service rate. 

Notice that the PDE in (|2.3|) with flow capacity constraint can be re-written as a scalar 
conservation law: 

d tP e (t, x) + d x c/ ) e (p e (t, xj) = (2.6) 

where the flux function satisfies (p e (p) = m.in{V e p, /x e |. 

A straightforward numerical scheme for solving system (12.3p - (12.5p is to apply space-time 
discretization and use the finite difference approximation [21 [T6| [2Tj 122] . For example, one 
could use an upwind scheme for conservation law (]2.3h and forward Euler method for ODE 
(|2.4p . moreover, to ensure numerical stability, one also needs to impose constraints on the 
time step, e.g. (jl.2p . 

One of the goals of this article is to provide an effective alternative to the numerical scheme 
mentioned above, namely the variational method or Lax formula [6l 1141 [T7] . The Lax formula 
was original proposed in [25] for the study of scalar conservation law of the form 

d t u + d x {f(u)) = (2.7) 

where the PDE (|2.3p is just a special case. The Lax formula expresses the entropy solution of 
(12. 7p as the solution of a minimization problem which, in our case, can be expressed explicitly. 
Notice that (12, 7h can be equivalently written as the Hamilton-Jacobi equation 

d t U + f(d x U) = 0, where U(t, x) = [ u{s,x)ds (2.8) 

Jo 



3 Variational method 

In this section, we apply the variational principle and reformulate the system governed by 
(I2.3p -( i2~3j) . we'll demonstrate the ability of the Lax-Hopf formula to handle both dynamics of 
the queue and processor, in a way consistent with ()2.3p - (|2.5p . and reduce the complexity of the 
coupling PDEs and ODEs introduced in [21J. To this end, we first consider a single processor, 
from now on, the superscript e will be dropped. We'll extend the following argument to general 
networks in Section 13.31 

Denote the flux Vp(t, x) = u(t, x), recall the density-based scalar conservation law (!2.6p . 
we rewrite it as a flow-based PDE (see [6] for more detail): 

d x u(t, x) + d t g(u(t, x)) = (3.9) 

where the function u i— > g(u) = p is defined as the inverse of the function p \— > <p(p) = u on 

the interval [0, p,/V], in other words, g(u) = mf{(/>(p) = u}, u € [0, p]. 

p 

Remark 3.1. The following analysis allows the relaxation ofg(-) to be any continuous, convex 
function. Similar technique can be applied to traffic models, see J<5|> 

The time horizon of consideration is [0, T] for any T > 0, and the processor is represented 
by a spacial interval [0, L\. Introducing the integral function {/(-, •) : [0, T] x [0, L] i— > R+ 

U (i, x) = / u(s, x) ds 
Jo 
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The conservation law (|3,9p can be equivalently written as a Hamilton-Jacobi equation, 

d x U(t, x) + g(d t U(t, x)) = (3.10) 

Define the cumulative goods that have arrived at the processor (possibly joining the queue) 
up to time t: Q(-) : [0, T] i-> E + . For what follows, Q(-) is assumed to be non-decreasing, 
continuous except for countably many t. We extend Q(-) to the interval (— oo, 0), with value 
assigned to it. For the rest of this article, we consider the left-continuous version of Q(-) 
where Q(t) = Q(t—). Given function Q(-), consider the Lipschitz continuous function: 

U(t) = inf (Q(r) + n(t- r)} < Q(t) (3.11) 

U (t) denotes the total goods that have been released from the queue to the processor up to 
time t. Notice that Q(t) — U(t) measures the length of the queue in front of the processor. 
We consider the following Cauchy problem (initial value problem by switching the roles of x 
and t): 

(d x ll(t, x)+g(d t U(t, x)) = 0, 
\U(t,Q) = U(t). 

For t > 0, the viscosity solution to ()3. 12[) is provided by the following Lax formula [6]. 
Proposition 3.2. The solution to 13.12\) satisfies 

"<■*•') = ^M^) +F «} = sMt) +(, I t ») (3J3) 

where g* is the Legendre transformation of g 

g*(u) = M{g(p)-up} (3.14) 

Proof. See 0. □ 
Proposition f)3.2|) has several impacts. (1) The first identity, known as the Lax formula, 
provides a grid free and accurate algorithm for solving the H-J equation. (2) The second 
identity suggests that the exit flow pattern can be obtained directly from the inflow profile 
Q(-), without intermediate computation of U(-), as a result, we no longer need two equations 
(|2.3p - ([2.4p to describe the queue and processor separately. (3) the assumptions on Q(-) suggest 
that the solution can accept very general boundary flows, for example, the inflow rate can 
exceed fi, or even be a distribution. 

3.1 Explicit instantiation of ( 13.13ft 

Due to the simple form of the flux function g(-) and its Legendre transformation <?*(■) consid- 
ered in this model, we're able to derive a explicit solution for the cumulative exit flow U(t, L) 
under minor assumption. To this end, we consider the family of piecewise affine functions Q(-), 
this assumption is not restrictive in application in the sense that all functions with certain 
regularity, e.g. piecewise continuous, can be approximated well by piecewise afflne functions. 
Notice that by (|3.1ip . U(-) is also piecewise affine. For a processor with length L, flow capacity 
H and processing speed V, the Legendre transformation g* in (|3.13p reduces to 



9*{u) 




U < y] 

U > y. 
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Setting x = L, formula (|3.13|) becomes 



U(t,L) = mini inf Q(r), inf \Q(t) - ut + (t - L/V)p,}\ . 

[r>t-L/V r<t-L/V J 

Observe that Q(-) is non decreasing and left-continuous, therefore 

inf {Q{j)-pr + {t-L/V)p] < lim (Q(r) - M r + (t - L/V)p) 

r<t-L/V T-+(t-L/V)~ 



(3.15) 



Q(t-L/V) = inf Q(r) 

r>t-L/y 



Now we can write (|3.15p in a reduced form 



U(t, L) = inf (Q(r) - /ur} + (t - L/V)n (3.16) 

r<t— L/V 

Remark 3.3. \3. 16\) has a simple interpretation, recall that the cumulative goods that have 
been released from the queue into the processor at time t is given by \3.11\) 

U(t) = M{Q(r) + p(t-r)} 

Then for a good that exits the processor at t, the time of its entry to the processor is t — L/V , 
thus the total goods that have exited the processor at time t is exactly equal to the total goods 
that have been released into the processor at time t — L/V , replacing t by t — L/V in 13.11\) 
gives \3.16\) . 

We also need to incorporate initial conditions into our formula, let q(t) be the queue length, 
p(t, x) be the density of goods on the processor, consider the following initial conditions 

q(0) = qo > 0, p(0, •) = po(0 €^([0, L]) (3.17) 
Let Q(-) : (— oo, T) i-> R + be the left-continuous function defined by 

x f° t < . 

Q(t) = < ~ (3.18 

\qo + Q(t) t > 

We're now ready to state and prove one of the main results in this article, namely, an 
explicit instantiation of the variational principle (|3.13p . with piecewise affine boundary profile 
Q(-) and initial conditions qo, Po{-)- 

Proposition 3.4. (Explicit formula) Given constants p., L, V and piecewise affine bound- 
ary datum Q(-) : [— oo, T) i— > M + with break points 6 1 : % € I}, we denote, for a given 
time t > L/V, the finite set Oj = {£j : £j < t — L/V}\^]{t — L/V}. Assume the following 
initial conditions 

9(0) = q , p(0, x) = p (x). (3.19) 

Moreover, define Q(-) as in \3.18\) , then the viscosity solution to VS. 12^ with initial condi- 
tions A3.19\) can be written as 

cL 

Po(()d( t G [0, L/V); 

U(t, L) = { ^-vt ^ 



L ~ vt (3.20) 



min (Q(r) - pr} + (t - L/V)p + / p (C) d( t G [L/V, T}. 

\t&h Jo 

U(-, L) is Lipschitz continuous. 
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Proof. Notice that L/V is the minimal time it takes to traverse the processor, therefore for 
< t < L/V, U(t, L) is completely determined by the initial distribution po(x) of goods 
inside the processor . A simple calculation using method of characteristics shows that for 
< t < L/V, 

U(t,L)= fv P {s,L)ds= [ t Vp (L-sV)ds = f p Q (() d( (3.21) 
Jo JO Jh-tV 

For t > L/V, we argue that a positive queue qo initially can be treated as an upward jump 
in Q(-) at t = 0, which is incorporated by replacing Q(-) with left-continuous, piecewise affine 
and non-decreasing function Q(-) defined in (|3.18p . Finally, recall (|3.16p . it remains to show 
that 

inf (Q(r) -fir} + (t-L/V)p = min (Q(r) - pr} + (t - L/V)p (3.22) 

T<t—L/V rEQt 

This is true since Q is piecewise affine, the infimum in (|3.22p must be obtained at either 
one of the break points, or at t — L/V, see Figure HJ Finally, since the Lax formula (13.131) 
gives viscosity solution to the Hamilton-Jacobi equation, U(t, L) is Lipschitz continuous with 
constant p. □ 




Figure 1: Graphical representation of the formula (|3. 16|) . The infimum of the vertical difference 
between two curves Q(j) and pr is obtained either at the break points of Q(-) (left), or at 
t - L/V (right). 

The solution H3.20n does not depend on any approximation and is therefore exact. Notice 
that the feasible set fit is finite, thus (|3.20p converts the variational principle from continu- 
ous form (I3.13P to discrete form assuming piecewise affine datum, this indicates a potential 
numerical improvement, see next subsection for more detail. 

3.2 Lax formula in discrete form 

Aim of this subsection is to derive the discrete version of (I3.20p for a given time grid. This 
new algorithm will be applied to network simulation and optimization later in this article. To 
this end, consider time horizon [0, T] for some T > and a uniform time grid = to < *i < 
. . . < ttf = T with step size h. 

Let Q(-) be any non-decreasing, left-continuous function, for shorter notation, let Qi = 
Q(U—), Ui = U(ti), < i < N. Then Q(-) is approximated by the piecewise- affine 
function with break points {(U, Qi)}j^ - To derive an explicit numerical scheme, we make one 
simplification by rounding ti — L/V to the nearest grid point to the left, and define integer 
constant A = \yj^~\ ■ A discrete version of Lax formula and the properties and error estimates 
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of the numerical solution is presented below. Note that the initial density distribution p$(x) 
solely induces integral terms in the Lax formula ()3.20p . whose error estimates are standard in 
current literature, and are of no particular interest to us. Therefore, in what follows, we only 
consider nonzero initial data qo. 

Proposition 3.5. (Lax formula in discrete form) Given constants L, V, \x, discrete values 
Qi, i = 0, . . . , N and initial data q(0) = qo > 0, po(x) = 0. Define Qo = 0, Qi = Qi+qo, 1 < 
i < N. Let U(t, L) satisfy i3.20\) . Then the following discrete version of Lax formula 

{0, < i < A; 

Q mm ^{Qj-fit^ + iti-L/V)^ A < i < N. (3 ' 23) 

satisfies 

< U hL < U i+1)L A < i < N-l (3.24) 

and 

< U i>L -U(U,L) < ([i-Ji-W A < % < N (3.25) 



Vh 1 Vh 

Vh 



In particular, if yj- is integer, i3. 23\) becomes exact. 



Proof. We first verify (13.241) . Notice that min iQj — utA— min (Qj — utA is equal 

0<i<i-A 1 J J ' 0<i<i+l-A 1 J J J 

to either or min (Qj — at A — (Qi+i-A — vU+i-a), in the latter case 
o<j<i-A ^ T T 

min {Qa — ut A — min {Qj — at A 
0<i<i-A L J J) 0<j<i+l-A L J J) 

< Qi-A - pti-A - (Qi+i-A - pU+i-a) < pA 

thus Ui l — Ui+i l = min {Qj — at A — min {Qj — atj] — Aa < 0. This shows 

o<i<«-A 1 J J ' 0<i<i+l-A 1 J J ' 

monotonicity. To show non negativity, it suffices to check that Ua l > 0. 

- L 
Vh 



For error estimate (|3.25p . recall A = \yr~\ , it follows from (13.20P and the definition of £lt 



that 



U(tj, L) = min (Q(r) - pr} + (U - L/V)p 

= min ljmn_ A {Q(tj) - ptj}, Q(U - L/V) - p(U - L/V) J + {U - L/V)p 

< u l>L 

and U(ti, L) <Uj L if and only if Q(tj - L/V) - p(tj - L/V) < min {Q{U) - fit A. In 

0<j<i— A 

this case, we have the estimate 

U i)L - Ufa, L) = 0< mm a {Q(tj) - ptj} - Q{U - L/V) + p{U - L/V) 

< Q{U-a) - pfi-A - Q(U - L/V) + n(u - L/V) 

< h(u-l/v-u-a) 

= n(U- L/V -ti + Ah) = (ih(A-^) 
This verifies (I335|) . □ 
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Remark 3.6. Proposition \3. 51 shows convergence of the numerical method in infinity norm 
\3.23\) as h — > 0, one can show similar result for more general Q(-) e.g. piecewise continuous. 

Since A > 1, the above proposition does not impose any constraint on the time grid, in 
terms of step size and uniformity. To improve the numerical precision of the approximation 
\3.23\) . by 13. 25\) . one may either increase the number of grid points or less obviously, choose 
h in a way such that the processing time L/V is a multiple of h. 

3.3 Extension to network 

Aim of this section is to extend the variational approach above to general supply chain network. 
Recalling Definition 12.11 for general network, we first introduce a few notations, for a fixed 
vertex v G V, the set of incoming arcs and outgoing arcs are denoted by T" and O v , respectively. 
In the case \O v \ > 1, we call the vertex v dispersive, and introduce the allocation rates A v ' e (t) 
for e G O" such that 

'0 < A v ' e {t) < 1, J2 A v,e (t) = 1, if|CF|>l 

eeo- (3.26) 
A v > e (t) = l, if|CF| = l 

These allocation rates A v > e (t) describe the percentages of flows coming from prior proces- 
sors going to the queue of downstream processor e, and are later subject to optimization. For 
a fixed processor e 6 i, let Q e (t) be a non decreasing, left-continuous function describing 
count of goods arriving at the buffer queue of the processor. Let W e {t) be a non decreasing 
Lipschitz continuous function describing the count of goods exiting the processor. Define con- 
stant processor length L e , processing speed V e and flow capacity . Adapting Proposition 
13.41 to network notations, we obtain a system of equations. 

Proposition 3.7. (Variational method for network model) Given a supply chain net- 
work, parameters L e , V e , y e , and initial conditions q e (0) = qo, p e (0, x) = Pq(x) for all e G A, 
then the following formulae provide solutions to this network model. 

j t Q e (t) = A v > e (t) ^ jW\t) a.e., v G V, e G O v (3.27) 
Q'iD ! ' 1 - ° (3.28) 




W e {t) 



-tv c 
inf 

{r<t-L e /V 



P e (()d(, te[0, L-/V e ) 

(Q e (r) - //Y} + (t - L e /V e )fi e + [ L p e (C) d(, t G [L e /F e , T] 

Jo 

(3.29) 

Remark 3.8. For any e G A, W e {t) is Lipschitz continuous by Proposition \ '3-4\ thus dif- 
ferentiable almost everywhere, therefore the right hand side of \3.21^ is well defined. As a 
consequence of {3.21^ , for e G O v , Q e is also Lipschitz continuous. 

The above system provides continuous-time explicit solutions to the network models orig- 
inally proposed in [21], where a system of coupling PDEs and ODEs were considered. As we 
shall see later, besides simplicity, (|3.29[) also yields better solution precision than the finite 
difference scheme for differential equations, under minor condition, making the variational 
method a competent candidate for large scale simulation and optimization. 
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4 Control and optimization on networks 



In this section, we look at a class of optimal control problems on networks: consider a general 
supply chain network given by Definition [21 for each of its processor and buffer queue, the 
dynamic is described by a Hamilton-Jacobi equation and the Lax formula. Assume at each 
dispersive node we can influence the allocation of the processed goods by a control A v ' e (t), in 
the goal of optimizing certain quantity of the network. A general cost functional is given by 

minVr(Q e ,r) (4.30) 

eeA 

The continuous-time optimal control problem is given by (j4.30|) . (I3.27|) . (|3.28|> and (I3.29p . 
It is shown later in this article that after a proper time discretization of the system and 
introducing binary variables, the discrete optimization problem can be formulated as a mixed 
integer program (MIP). We'll also discuss several model extensions afterwards. 

4.1 Time discretization 

We adapt the "discretize-then-optimize" strategy by first time-discretize system (|3.27j) - (|3.29[) 
and then optimize the discrete system. To do this, assume a uniform time grid {ti}^L with 
step size h, and use notation Q\ = Q e (ti), Wf = W e (ti), Af e = A v ' e (ti), A e = \J4 K \. To 
avoid differentiation in (|3.27p . we use the following lemma. 



Lemma 4.1. Equation (3.21) can be equivalently written as 

£ Q e (t) = y, w ~ e V) ( 4 - 31 ) 

Q e (h) < Q e (t 2 ), for any h < t 2 , e G 0\ (4.32) 
Proof. On one hand, by (|3.26p . summing (|3.27|) over e £ O v gives (|4.31|) . Moreover, 

j t Q e (t) = A">°(t) Y, jW\t) > a.e. 
this establishes (|4.32p . On the other hand, differentiating (|4.31|) gives 

define A v ' e (t) = f t Q e {t) / 'Y^sai W%t), e £ O v , then it's straightforward to verify that t$SM) 
holds. □ 
Note that (|4.3ip (|4.32p guarantee conservation of flux and non negativity of flow. After 
applying Lemma [4. 11 information of the allocation rates A v ' e is hidden in Q e and it is a matter 
of simple calculation to recover them from the solution. The discrete version of (|3.28p . (|3.29p 
becomes 

Qi = Qo + Qt, 1 < i < N (4.33) 

i < A e 

W * = t min {Qf -//%} + (U - L e /V e )n e i > A e (434) 

0<j<«-A e 
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To overcome the difficulty originated from the nonlinear function "min" in (|4.34|) , we write 
(|4.34p in a slightly different way. Introducing variables R\ for e £ A, A e < i < N, defined 
recursively by 

Re _ fQi_A«-/* 6 *i-A«, » = A e 

* jminji^, Q?_ A e-M e ti-A«}, A e <i<N + A e [ ' ' 

Then the second equation of (I4.34p becomes 

Wt = R\ + (U- L e /V e )n e , i > A e (4.36) 

By virtue of this new variable Rf, for each time step i, we need to evaluate the "min" function 
only once, this reduces the number of operations and variables in the problem. The next 
step leads to the mixed integer formulation: introduce binary variables /3f 6 {0, 1}, i = 
A e + 1, . . . , A e + N, condition (]4.35p is equivalent to 



R%e — Qq — //to; ' 



{ RU + - 1)M < Rf < RU 

(4.37) 

LQ?_ Ae -)U%-A«-M^ < Rf < Q?_ Ae - fl%. A e 



where constant M is chosen to be a sufficiently large number. With all preliminary results at 
hand, we're now ready to state the main result of this article. 

Theorem 4.2. (MIP formulation of optimal control problem) 

Given a supply chain network as in Definition \2.1l and parameters L e , V e , fi e , e € A, 
consider initial conditions q e (0) = q$, e € A and quantities measured by 
(j 7e (Q e , W e , -R e )) eg ^. Apply a uniform time discretization = to < t\ . . . < i/v = T, with step 
size h, then the optimal control problem \4--30 ), (3.21), \3. 28}) and \3. 29\) can be formulated 
as a mixed integer program. 

min V T e (Q e , W e , R e ) (4.38) 

subject to 

E Qi = E W ^ Qi-i - ®i> l<^<N (4.39) 

R e Ae = -fj, e t , R\_ x + (/3 4 e - 1)M < Rf < R\_ x , A e + l<i<A e + iV (4.40) 
Qo + Q!-Ae ~ H%-A* - Mffi < R\ <q e + QI_ Ae - fJ,%- A «, A e + l<i<A e + N (4.41) 
Wf = R\ + (U - L e /V e )n e A e < i < N (4.42) 

A e = r^l, Pfe{0, 1}, e e A, vGV. (4.43) 
where M is a sufficiently large constant. 

Proof. We have already established (|4.39|) and (I4.42p in Lemma 14.11 and f|4.36j) , respectively. 

(|4.40p . (|4.4ip are immediate consequences of (|4.33p and (|4.37p . □ 
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Remark 4.3. If \O v \ = 1, i.e. node v is non-dispersive, the non-negative flow condition 
Qf-i — Qf i n 1^4-39 ) is automatically satisfied since Wf, e£l" are guaranteed to be non- 
decreasing w.r.t. time step i, by Proposition \3. 51 



4.2 Model extensions 

1. Finite size buffers. In the design of production lines, it is mandatory to limit the size 
of the buffering queues. To implement this constraints in the MIP formulation, we need to 
express the queue lengths explicitly. Recall (|3. 11 [) that U (•) describes the count of goods that 
have been released from the queue into the processor. The queue length at time t is clearly 

q e (t) = Qe(t)-U e (t) = (Q e (t)-^t)-inf{Q e (r)- M e r} (4.44) 

apply time discretization i = 0, . . . , JV, 

qf = Ql-ifti - min.{Q|-/i%} = Q\ - - Rf +Ae (4.45) 

0<j<i 

The finite size buffer condition can then be implemented in the mixed integer formulation by 
adding the following linear inequalities 

Qf - fx% - R e i+A e < const, e G A, 0<i<N. (4.46) 

2. Optimal inflow profile. It is sometimes necessary to consider costs for storing and 
managing surplus parts waiting in the queue. The problem now arises as how to choose the 
inflow profile to the network such that the cost from queuing is minimized while a satisfactory 
production of goods is guaranteed. This can be achieved by adding the following term to the 
cost function (14.381) 

JV 

E °\ YM " - Rf+A°), const c\ > 0, e G A (4.47) 

eeA i=0 



4.3 Comparison with existing approaches 

The MIP (|4.38p - (|4.43p differs from the one in |20| in several aspects. First, the governing 
equations for the system in this paper is derived from a variational principle, the main vari- 
ables are cumulative good counts, which are connected via a Hamilton-Jacobi equation. The 
Lax formula governs dynamics of both the queue and the processor, and avoids issues like 
discontinuous dependence on the queue induced by (|2.5p . and in some rare cases, non-physical 
solutions (e.g., see Section [575]) . Our approach accepts boundary datum Q{-) that are dis- 
continuous, which implies a Dirac delta in the flow, such situation cannot be handled by the 
conservation law without a separate modeling of the queue, this intuitively explains why in 
our formulation, an explicit treatment of the queue is unnecessary. 

Second, the Lax formula provides closed-form and exact solution to the system under 
minor conditions, it does not rely on any spacial variable or finite difference approximation, 
this reduces memory usage and the number of flops compared to existing numerical schemes 
[20] l2l] . Moreover, the approach proposed here does not impose constraints on the time step 
size for the sake of numerical stability; on the other hand, the finite difference scheme usually 
yields larger system due to CFL condition for the PDE and stiffness condition for the ODE 
[20]. 
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Next we compare the MIP (|4.38|) - (|4.43|) with the one presented in |2Uj with a two-point 
coarse spacial discretization, in terms of number of variables and equality/inequality con- 
straints. Assume the numbers of arcs and nodes in the network are \A\ and |V| respectively, 
and the number of time intervals is N. Then the MIP proposed in [20J consists of 3iV|^4| (real) 
+N \A\ (binary) variables and (2|«4| + |V|)iV (equality) +10|„4|iV (inequality) constraints. On 
the other hand, the MIP in Theorem 14.21 has 3iV|„4| (real) +N \A\ (binary) variables, and 
(|-4| + |V|) N (equality) +7\A\ N (inequality) constraints. 

Although the number of (binary) variables and constraints in both MIP formulations seem 
comparable given the same network and time grid, in order to obtain certain accuracy of the 
solution, a two-point upwind discretization of the PDE ([20]) is insufficient, therefore finer 
spacial grid and, as a result of CFL condition, finer time grid have to be used; on the other 
hand, due to the good precision of the variational method, we are able to achieve the same 
level of accuracy by using much less time grid points and no spacial grid points. As we shall 
see later in Section 15.41 (|4.38p - (|4.43l) can yield MIP that is significantly smaller in size than 
the MIP proposed in [20]. 

5 Numerical example 

In this section, we conduct a sequence of numerical studies of the variational approach and 
the resulting MIP, we also incorporate model extensions discussed in Section 14.21 Consider a 
network example in Figured! it consists of seven arcs (processors) a — g and four nodes 1 — 4. 
There are two control variables A l,b (t) and A 2,e (t) at dispersive nodes 1 and 2, respectively. 
Set the time horizon [0, 10], and assume an inflow f a {t) into processor a given by f a {t) = 
37.5, < t < 2, f a (t) = 0, otherwise. Network parameters are shown in Table 1. In addition, 
we assume zero initial conditions q$ = 0, Pq(x) = 0, e € V. 




Figure 2: A network example consisting of seven arcs and four nodes 



Processor 
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c 
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/ 


9 


L e 


2 


2 


2 


2 


2 


2 


2 


V e 


2 


1 


2 


4 


2 


2 


2 


P e 


15 


6 


5 


4 


3.5 


8 


14 



Table 1: Processor parameters 



All mixed integer programs presented in this section are solved with ILOG Cplex 12.1.0 
[24J, on the HPC platform with Intel Xeon X5675 Six-Core 3.06 GHz processor, provided by 
the Penn State High Performance Computing Systems [31] . 
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5.1 Without buffer size constraints 



Fix the inflow f a (t), and assume infinite buffers sizes, the objective is to maximize the total 
throughput of the network within the time horizon, namely W 9 (10). This is achieved by 
using the MIP in Theorem 14.21 , with a time grid of 200 points. The optimal allocation 
rates A 1,b (t), A 2,e (t) are illustrated in Figure [3] and HI In the optimal network solution, only 
processor a, b and c possesses non-empty queue, as shown in Figured! The optimal throughput 
W 9 (10) is 58.75. 
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Figure 3: Without buffer size constraints: opti- Figure 4: Without buffer size constraints: opti- 
mal allocation rate A 1 ' b mal allocation rate A 2 ' e 
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Figure 5: Without buffer size constraints: queues in front of each processor 



5.2 With finite buffer constraints 

Consider a similar optimal control problem as Section I5.ll except for an upper bound for the 
queue sizes, namely, q% < 10, q\ < 10, i = 0, . . . , N. This is implemented by adding inequality 
constraints (|4.46p to the MIP. The queues in the new optimal solution are shown in Figure [HJ 
Notice that the queue in front of processor a cannot be controlled once the inflow f a is fixed. 
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Figure 6: With buffer size constraints: optimal Figure 7: With buffer size constraints: optimal 
allocation rate A 1 ' b allocation rate A 2, e 
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Figure 8: With buffer size constraints: queue in front of each processor. 

With the constraints on the queues, the adjustment of the allocation rate A^' b compared 
to the infinite buffer case can be clearly depicted in Figure EJ see also Figure [3] for compari- 
son. Interestingly, the optimal value of VF S (10) is 58.75, same as the one without any queue 
constraints. This implies the non uniqueness of the optimal solution. 



5.3 Optimal inflow profile 

In view of model extension in Section 14.2} we consider extra cost for queuing by applying the 
following objective function, 



N 

max e {^(10) - £ c\ £(Q? - i?f +Ae )} 



Of, A 



(5.48) 



where the unit queuing cost c e q = 1 for all e £ A. Treating the inflow profile Q a as free variable, 
we want to maximize the production of the network while keep network congestion (queue) to 
the minimum. This is implemented in the MIP (14.38p -( l4.43|) with cost (I5.48p . The resulting 
optimal inflow is shown in Figure [91 optimal allocation rates at vertices 1 and 2 are illustrated 
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in Figure [TU1 and [TT1 Under the optimal inflow and these allocation rates, the queues in front 
of every processor are empty throughout the time horizon. The final throughput W 9 (10) is 
again 58.75. 
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Figure 9: Optimal inflow into arc a. 
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Figure 10: Optimal inflow case: allocation rate Figure 11: Optimal inflow case: allocation rate 

A l,b A 2,e 



5.4 Solution Precision 

In the remainder of Section [5] we evaluate the new algorithm and the resulting MIP in terms of 
solution precision and run time, and compare with the finite difference scheme (TBI 1201 12T1 [22] 
and the corresponding MIP [20]. To confirm the assertion we made at the end of Section^ i.e. 
with the same level of solution precision, the MIP (j4.38[ )- (14.43|) can be significantly smaller 
in size than the MIP in [20J, we use the same network and parameters as previous examples, 
extend the time horizon to [0, 80] and assume infinite capacity of queue buffers. The objective 

is again to maximize the total throughput max Wi T . The inflow is fixed 

Qf,W? N 

fit) = < - ° - 1 " 10 (5.49) 

] - 10 < t < 80 
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It is expected that the total throughput W 9 (80) is equal to 450 provided the time hori- 
zon is large enough. We run instances of the mixed integer program (|4.38p - (|4.43p for N = 
40, 50, 60, . . . , 1000 where N is the number of time intervals. The optimal throughputs W% 
are plotted in Figure [T2j 




Figure 12: Optimal output for different choices of time grid 

We interpret Figure [12] as follow: recall Proposition 13.51 that the discrete Lax formula is 
exact provided that the throughput times L e /V e , e S A are multiples of time step h, such 
as cases N = 160, 320, 480, 640, 800, 960 in this numerical example. Otherwise, the error is 
given by estimate (|3.25 j> . 

o < w K -wm < (r^i-^> 

thus the numerical values > the exact value W 9 (80) = 450, and displays an oscillatory pattern 
with damping, due to the discontinuous function \yr \ — and the factor h. Notice that the 
right hand side of error estimate (|3.25|) can be further relaxed to hji, one can easily prove the 
upper error bound for network throughput: 

< W 9 N -W 9 (T) < /(max^/i e (5.50) 

P ^ eGp 

where p = {ei, e2, . . . , e m } is any viable path, ej G A, W 9 (T) is the continuous-time solution. 
In view of the network and Table 1, the upper bound for the error is plotted in Figure [T2l 

We also implement the MIP [20J based on a finite difference approximation of the system, 
under the same numerical setting. We choose different values of time intervals ./V and spacial 
intervals D, for example, D = 1 means the two-point discretization. The results are also 
plotted in Figure PT2l notice when D increases, the minimal value of N allowed by CFL 
condition also increases. This numerical experiment leads to two observations: 1. the values 
of D £ {1, 2, 4} seem to have no effect on the error as long as the CFL condition is satisfied, i.e. 
N large enough; 2. in terms of solution precision, the worse case in our approach outperforms 
the finite difference approach. 

The above numerical experiment suggests that the size of the MIP can be significantly 
reduced using our formulation, while attaining certain solution precision. For example, to 
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control the variation of the total throughput of the network within [450, 452], a time grid 
of over 2000 points is required for finite difference-based approach, while for the variational 
approach, a time grid of N = 80 or 160 is sufficient to achieve such solution precision. 



5.5 A case study of a smoothed out version of (12.51) 

In the formulation of [20j, in order to simplify the ODE and avoid numerical difficulty such 
as discontinuity in the state variable, a smoothed out version of the ODE (I2.4p - (|2.5p was used 
instead 

d : q\t) = u e (t)-r{p e (t,a e )) (5.51) 



dt 

f e {p e (t, a e )) 



mm < [i 



'it) 



(5.52) 



where e is a smoothing parameter. This approximation has its merit in many computational 
practices, there are however, a few cases where this could lead to unwanted solutions. For 
instance, let the network inflow of previous example in Section 15.41 be 



nt) 



< t < 10 
10 < t < 80 



(5.53) 



the MIP in [20J is solved for N = 600, D = 1. The inflow, exit flow and queue length of 
processor a are plotted in Figure [T3l A forward discretization of (|5.5ip - (|5.52p yields a queue 
of non-negligible size, even if the inflow is strictly below the processor capacity fi a = 15. We 
argue that such unwanted queue may be reduced by choosing smaller e. But the numerical 
stiffness condition ([20]) requires At < e, this is again a trade-off between problem size and 
solution quality. On the other hand, the variational method handles these problems well with 
N = 600, see Figure M 





Figure 13: Finite difference with e = 0.5 



Figure 14: Variational approach 



5.6 Solution time 

As our final experiment, we test the numerical efficiency of two MIP approaches using the 
same network without buffer size constraints. The time horizon is set to be [0, 10], fix the 
network inflow to be 




< t < 5 
5 < t < 10 
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For the MIP presented in this paper, the objective function is chosen to be maxW^; for the 
MIP in [20j, we use a two-point spacial discretization, i.e. D = 1, the objective function is 
adjusted to niax^^ 1 Wi/(1 + ij), where u>i is the exit flow on processor g at time t$. 

We run the two MIPs for number of time intervals N ranging from 160 to 3000, the 
recorded solution times are summarized in Figure [T5l Although the two MIPs are similar in 
size, the run times of our MIP problem is significantly lower than its opponent, for N > 1500, 
and we observe a nearly linear growth of solution time vs. N. Possible mechanism leading to 
such significant difference for large scale MIP problems, either involving internal structure of 
the system discretization or specific settings of the branch-and-bound algorithm, is currently 
under investigation. 




Figure 15: MIP computational time 



6 Conclusion 

In this article we propose a variational method for modeling and formulating continuous supply 
chain networks. The solution is derived in continuous/discrete time and the formulation is 
independent of the partial differential equations which is arguably more general than the 
approach in [20, [21] in terms of modeling assumption and solution method. We demonstrate 
the unique features of the Lax- formula such as explicitness, grid- free and exactness and verify 
them numerically. 

We also apply the variational approach to network optimization and end up with mixed 
integer programming. By appropriately choosing the time grid, the MIP proposed in this 
article is much smaller in size and can achieve better solution precision than the MIP in [20!. 
In a comparison of numerical efficiency , we also observe, in certain case, significant speed-up 
by using our MIP approach. 

Future research includes a continuous formulation of the optimality condition for system 
(|4.30p . (|3.27p . (|3.28p and (|3.29p . which is viewed as a bi-level optimization problem; and the 
incorporation of more realistic conditions such as uncertainty, nonlinearity and inhomogeneity 
of the production process. 
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